R como una plataforma SIG

El análisis espacial es una extensión natural a las capacidades de R. Existen múltiples paquetes para manipular y analizar datos espacio-temporales en R. En esta clase vamos a dar una introducción a esto.

A la hora de llevar a cabo visualizaciones y análisis espaciales en R se heredan todos los beneficios que se han venido mencionando en clases sesiones pasadas.

El contra más grande de usar R pos sí solo como SIG

Por esto es recomendable parear R con software SIG especializado, por ejemplo QGIS, que es sumamente poderoso y de código abierto:

http://www.qgis.org/en/site/

Formatos de datos espaciales

En el mundo del análisis espacial existen dos formatos principales:

El formato vectorial, descripción y ejemplos

El formato vectorial se refiere a estructuras que se conforman por puntos, líneas o polígonos asociados a una tabla de datos. Por ejemplo poner unos ejemplos:

Cada instancia de un crimen asociado a un lugar (un par de coordenadas) se puede asociar a un formato vectorial de puntos. Por ejemplo un mapa de Homicidios dolosos en la CDMX desde Enero del 2013 se podría visualizar de la siguiente manera usando un formato vectorial:

Las carreteras y caminos son candidatos naturales para digitalizarse como líneas.

cada delimitación política de un estado de la república es un polígono, por lo tanto una manera natural de almacenar digitalmente un mapa de los estados es a través de un archivo vectorial. Como se mencionó, estas unidades (polígonos en este caso) pueden estar asociadas a una tabla de datos, por ejemplo, las variables extensión total en Ha y población total.

Un primer objeto espacial en R e introducción a APIs

Casi todas las estructuras de datos vectoriales en R se basan en el paquete. Hoy en día es prácticamente un estándar y la mayoría de los otros paquetes construyen sobre este.

El paquete sp está construído de manera tal que está dividido en funcionalidades que permiten crear y manipular los objetos de la sección anterior: puntos, líneas y polígonos.

Para incentivar el entendimiento de estas estructuras de datos, construyamos un objeto espacial de cero usando el paquete sp.

El geocoding es el proceso computacional de transformar una dirección (texto) a coordenadas geográficas. Existen numerosos sevicios para llevar esto a cabo, el más popular es el de Google Geocoding API:

https://developers.google.com/maps/documentation/geocoding/start

API son las siglas en inglés para Application Programming Interface. Una API no es lo mismo que un servidor remoto, es una parte de un servidor que recibe peticiones y regresa respuestas, esto es; una API es una manera de ofrecer servicios a clientes. En este caso, Google, a través de la API anterior ofrece el servicio de geocoding.

Para usarla se debe generar y activar una key en esta dirección:

https://developers.google.com/maps/documentation/geocoding/get-api-key

Luego conviene almacenarla como un objeto en R.

key_google_geocoding = "AIzaSyCup9GDTLlNdZgc9F_sUZjObp4UynwmNNM"

En R existen numerosos paquetes para llevar a cabo consultas de geocoding, a continuación usaremos el paquete “googleway” para obtener datos (entre ellos las coordenadas) de un conjunto de cines cinépolis.

library("googleway")
## Warning: package 'googleway' was built under R version 3.4.2
resultados_cinepolis = google_places(search_string = "cinepolis,CDMX,Mexico",
                                place_type="movie_theater",key = key_google_geocoding) # se debe de utilizar la key anterior

names(resultados_cinepolis)
## [1] "html_attributions" "next_page_token"   "results"          
## [4] "status"

La función básica “names()” nos regresa los nombres de los elementos de una lista o de las columnas de un data.frame. Usándola sobre el objeto anterior, vimos que en la lista se encuentran los objetos “html_attributions”, “next_page_token”, “results” y “status”.

El objeto “results” es un data.frame que tiene todos los datos sobre las direcciones encontradas al llevar a cabo el proceso de geocoding. De estas vamos a elegir un subconjunto útil.

library("tidyverse")
## Loading tidyverse: ggplot2
## Loading tidyverse: tibble
## Loading tidyverse: tidyr
## Loading tidyverse: readr
## Loading tidyverse: purrr
## Loading tidyverse: dplyr
## Warning: package 'dplyr' was built under R version 3.4.2
## Conflicts with tidy packages ----------------------------------------------
## filter(): dplyr, stats
## lag():    dplyr, stats
cinepolis = data.frame(select(resultados_cinepolis$results,place_id,name,formatted_address,rating),
                       lat = resultados_cinepolis$results$geometry$location$lat,
                       lon = resultados_cinepolis$results$geometry$location$lng,
                       stringsAsFactors = FALSE)

El data.frame “cinepolis” todavía no es un objeto espacial. Aún así, utilizando la función “gvisMap()” en el paquete “googleVis” podemos visualizar estos cines de manera dinámica. El parámetro esencial de esta función es el nombre de la columna donde se encuentran las coordenadas de los puntos, en esta columna las coordenadas deben estar en formato texto y con la forma:

lat:lon

construyamos un objeto con esta forma y agreguémoslo como columna a nuestro data.frame cinepolis.

cinepolis$google_latlon = paste0(cinepolis$lat,
                                 rep(":",nrow(cinepolis)),
                                 cinepolis$lon)

Ya que este data.frame tiene esta columna en el formato apropiado podemos utilizar la función “gvisMAP()”, el argumento “tipvar” indica qué info estará asociada a los marcadores en el mapa que desplegaremos.

library("googleVis")
## Warning: package 'googleVis' was built under R version 3.4.2
## Creating a generic function for 'toJSON' from package 'jsonlite' in package 'googleVis'
## 
## Welcome to googleVis version 0.6.2
## 
## Please read Google's Terms of Use
## before you start using the package:
## https://developers.google.com/terms/
## 
## Note, the plot method of googleVis will by default use
## the standard browser to display its output.
## 
## See the googleVis package vignettes for more details,
## or visit http://github.com/mages/googleVis.
## 
## To suppress this message use:
## suppressPackageStartupMessages(library(googleVis))
cinepolis_en_google <- gvisMap(cinepolis, locationvar="google_latlon", tipvar="name", 
              options=list(showTip=TRUE, showLine=F, enableScrollWheel=TRUE, 
                           mapType='satellite', useMapTypeControl=TRUE, width=600,height=400))

# luego las funciones 

# plot(cinepolie_en_google) abre una ventana con el mapa

# print(cinepolis_en_google, file="tu_ruta/cinepolis.html") guarda el mapa en un archivo html

¿Qué es una proyección cartográfica? - objetos espaciales en R

Ahora usaremos los paquetes “rgdal” y “sp” para convertir el data.frame anterior a uno espacial usando la función “coordinates()” para indicarle a R qué columnas serán sus coordenadas. Esto convierte al objeto “cinepolis” en un data.frame espacial.

library("sp")
library("rgdal")
## Warning: package 'rgdal' was built under R version 3.4.2
## rgdal: version: 1.2-15, (SVN revision 691)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.0, released 2017/04/28
##  Path to GDAL shared files: D:/software/R/R-3.4.1/library/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: D:/software/R/R-3.4.1/library/rgdal/proj
##  Linking to sp version: 1.2-5
coordinates(cinepolis) =~ lon + lat

class(cinepolis)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

Proyectar es llevar a cabo un proceso para hacer que un objeto sea visible sobre otro. En geografía, una proyección se basa en un modelo de la tierra para mapear coordenadas de la misma hacia un plano (x,y).

En este momento nuestra data.frame cinepolis no está proyectada:

summary(cinepolis)
## Object of class SpatialPointsDataFrame
## Coordinates:
##           min       max
## lon -99.25832 -99.01790
## lat  19.29530  19.50611
## Is projected: NA 
## proj4string : [NA]
## Number of points: 20
## Data attributes:
##    place_id             name           formatted_address      rating     
##  Length:20          Length:20          Length:20          Min.   :3.700  
##  Class :character   Class :character   Class :character   1st Qu.:4.075  
##  Mode  :character   Mode  :character   Mode  :character   Median :4.200  
##                                                           Mean   :4.220  
##                                                           3rd Qu.:4.400  
##                                                           Max.   :4.600  
##  google_latlon     
##  Length:20         
##  Class :character  
##  Mode  :character  
##                    
##                    
## 

Existen muchas proyecciones, la oficial para nuestro país (INEGI) es una basada en un modelo cónico de la tierra. Esta refiere elementos que se hallan en la esfera de la Tierra sobre un cono tangente, usando de vértice el eje que une a los dos polos.

Una proyección particularmente común es WGS84 latitude-longitude, por ejemplo es la utilizada por google maps.

un Coordinate Reference System (CRS) es la combinación de un sistema de coordenadas geográficas y posiblemente de una proyección. Ya que tenemos nuestros coordenadas asignemosle una proyección. Esto se puede hacer de múltiples maneras, las proyecciónes populares tienen asignado un código:

  • WGS84 latitude-longitude := CRS(“+init=EPSG:4326”)

O se puede definir parámetro por parámetro:

  • WGS84 latitude-longitude := CRS(“+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs”)

Algunos de los parámetros son evidentes, por ejemplo +ellps=WGS84 se refiere a que el elipsoide 3d que se usa para modelar a la tierra es uno llamado WGS84. Sabemos que las coordenadas serán longitud y latitud. Aún así, entender los parámetros de distintas proyecciones escapa el alcance de esta clase.

Lo más importante aquí es saber que si se quieren asociar espacialmente dos conjuntos de información (por ejemplo dos mapas vectoriales), deben encontrarse en la misma proyección.

Asignemos a nuestro data.frame la proyección WGS84 latitude-longitude, esto se puede hacer con la función proj4string()

# asignamos una proyección a estos puntos

proj4string(cinepolis) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

Ahora, este objeto está bien preparado para escribirse en disco y usarse como una verdadera capa SIG.

#writeOGR(cinepolis, "cinepolis.shp", "cinepolis", driver="ESRI Shapefile")

Este archivo se puede cargar nuevamente como un data.frame espacial con la función “readOGR()” en el paquete “rgdal”

# ruta_shapefile = "mi_ruta/cinepolis.shp"

# data_frame_espacial = readOGR(ruta_shapefile,ogrListLayers(ruta_shapefile)[1])

Cargar el shapefile generado y con base en él volver a generar el mapa con la función “gVisMap()” nota que esta función recibe un data.frame normal, no uno espacial.

Un data.frame espacial está compuesto por una tabla de datos asociada a una estructura espacial (geometría). Se puede acceder a la tabla de datos haciendo

# data.frame@data

O se puede convertir el data.frame espacial a uno usual con la función “as.data.frame()”.

Por otro lado hay que tener cuidado con los nombres de las columnas puesto en el proceso de escribir un shapefile a disco generalmente se acortan los nombres de columna muy largos.

Cómo usar el DENUE y no morir en el intento

Las funciones de geocoding no están diseñadas para generar datos exhaustivos sobre un negocio en particular. De hecho, el límite de resultados que te puede arrojar una petición en particular es de 20, esto 3 veces utilizando el page_token para un total de 60 resultados (ver la ayuda de la función “google_places()”).

Existen trucos para sobrellevar esta limitante pero en general no es la mejor manera de abordar el problema de ubicar exhaustivamente negocios, por ejemplo haciendo algo como lo que se presenta aquí (aunque usan el lenguaje python se puede implementar en R):

https://iliauk.com/2015/12/18/data-mining-google-places-cafe-nero-example/

En México existe una base de datos decente de negocios. El DENUE, que ahora está asociado a una API y al paquete de R “inegiR” para hacer consultas desarrollados por INEGI mismo.

Para utilizar este paquete de R, análogamente a como se hizo con la API de Google Geocoding se debe generar una key (nota que el DENUE tiene su propio generador de keys):

http://www.beta.inegi.org.mx/servicios/api_biinegi.html

Por supuesto también es buena idea guardarla en un objeto de R.

# asignar key a un objeto
key_denue = "b93b9e6c-5ec5-482d-bce4-03321c402f70"

Hay dos maneras de trabajar con el paqueteR para obtener datos del INEGI:

http://www.beta.inegi.org.mx/servicios/api_biinegi.html

Por suerte, los datos del DENUE se encuentran en este segundo caso. Por ejemplo la función “denue_inegi()” permite hacer una búsqueda espacial alrededor de un punto de referencia. Por ejemplo alrededor de un cine de nuestra elección.

library("inegiR")
## Warning: package 'inegiR' was built under R version 3.4.2
# elijamos un cine
cine = as.data.frame(cinepolis[2,])

# luego hacemos la consulta espacial buscando por ejemplo cines cinemex
cinemex_alrededor <- denue_inegi(latitud = cine$lat, 
                                  longitud = cine$lon, 
                                  token = key_denue,
                                  metros=5000,
                                  keyword="cinemex")

cinemex_alrededor
##        id                Nombre                         Razon
## 1 6309618      CINEMEX TEZONTLE OPERADORA DE CINEMAS SA DE CV
## 2 6314053 CINEMEX PLAZA ORIENTE OPERADORA DE CINEMAS SA DE CV
##                                                  Actividad
## 1 Exhibición de películas y otros materiales audiovisuales
## 2 Exhibición de películas y otros materiales audiovisuales
##            Estrato Vialidad             Calle NumExterior NumInterior
## 1 31 a 50 personas  AVENIDA CANAL DE TEZONTLE        1512            
## 2 31 a 50 personas  AVENIDA CANAL DE TEZONTLE        1520            
##                   Colonia   CP                                Ubicacion
## 1  DR ALONSO ORTIZ TIRADO 9020 IZTAPALAPA, Iztapalapa, CIUDAD DE MÉXICO
## 2 DR ALFONSO ORTIZ TIRADO 9020 IZTAPALAPA, Iztapalapa, CIUDAD DE MÉXICO
##   Tel eMail SitioWeb Tipo     Longitud     Latitud
## 1                    Fijo -99.08240723 19.38199430
## 2                    Fijo -99.08169257 19.38182559
# podríamos estar interesados también en otros cines menos conocidos

# para estos se puede primero buscar negocios cuyo nombre incluya 
# la palabra cine
negocios_alrededor <- denue_inegi(latitud = cine$lat, 
                                  longitud = cine$lon, 
                                  token = key_denue,
                                  metros=5000,
                                  keyword="cine")

# luego descartamos todos los negocios que no sean cines
cines_alrededor = filter(negocios_alrededor,
                         Actividad=="Exhibición de películas y otros materiales audiovisuales")

cines_alrededor
##        id              Nombre                           Razon
## 1 1031205       CINE NACIONAL                                
## 2  891138 LA CASA DEL CINE MX CONSULADO CINEMATOGRÁFICO, A.C.
## 3 1001689          CINE VENUS         DISTRITUR, S.A. DE C.V.
##                                                  Actividad         Estrato
## 1 Exhibición de películas y otros materiales audiovisuales 6 a 10 personas
## 2 Exhibición de películas y otros materiales audiovisuales 6 a 10 personas
## 3 Exhibición de películas y otros materiales audiovisuales  0 a 5 personas
##   Vialidad                        Calle NumExterior NumInterior Colonia
## 1  AVENIDA FRAY SERVANDO TERESA DE MIER           0              CENTRO
## 2    CALLE         REPÚBLICA DE URUGUAY          52           2  CENTRO
## 3    CALLE           REPÚBLICA DE CHILE          34              CENTRO
##      CP                                Ubicacion        Tel
## 1 06090 CUAUHTÉMOC, Cuauhtémoc, CIUDAD DE MÉXICO 5555224024
## 2 06000 CUAUHTÉMOC, Cuauhtémoc, CIUDAD DE MÉXICO 5555124243
## 3 06020 CUAUHTÉMOC, Cuauhtémoc, CIUDAD DE MÉXICO 5555231099
##                              eMail             SitioWeb Tipo     Longitud
## 1 ESTRELLA.MARINLOPEZ@YAHOO.COM.MX                      Fijo -99.12994101
## 2        CONTACTO@LACASADELCINE.MX WWW.LACASADELCINE.MX Fijo -99.13762812
## 3 ESTRELLA.MARINLOPEZ@YAHOO.COM.MX                      Fijo -99.13563358
##       Latitud
## 1 19.42352399
## 2 19.43067352
## 3 19.43788506

Por supuesto lo anterior nos permite automatizar el proceso de búsqueda de otros cines alrededor de cines cinepolis.

Una breve introducción a iteración

Una herramienta útil para evitar escribir más código son los mecanismos de iteración. Estos son útiles cuando se quiere llevar a cabo la misma tarea múltiples veces.

En esta clase aprenderemos sobre un paradigma de iteración: la programación imperativa.

Este paradigma es el más antiguo y permite introducirse fácilmente al tema pues hace que la iteración sea muy explícita. Como desventaja, las estructuras de este paradigma, llamados bucles (en inglés loops), tienden a ser verbosos.

Existe otro paradigma llamado programación funcional que ofrece herramientas para extraer todo el código duplicado para que cada bucle tenga su propia función. Después de un poco de práctica se puede resolver los problemas más comunes en iteración de manera más sencilla, utilizando menos código y por lo mismo cometiendo menos errores. Revisar el paquete “purrr”:

https://cran.r-project.org/web/packages/purrr/index.html

Programación imperativa

Un bucle está compuesto de una secuencia y un cuerpo que se ejecuta con base en esta secuencia.

Una secuencia:

1:5
## [1] 1 2 3 4 5

Un bucle for puede ejecutar el mismo código (el contenido en su cuerpo) variando un iterador. Por ejemplo el siguiente código nos muestra uno por uno los números en la secuencia misma y en cada una de esas iteraciones nos muestra la palabra “gatito”.

for (i in 1:5) # para la secuancia del 1 al 10
{
  print(i)
  print("gatito") # ejecutar esto 
}
## [1] 1
## [1] "gatito"
## [1] 2
## [1] "gatito"
## [1] 3
## [1] "gatito"
## [1] 4
## [1] "gatito"
## [1] 5
## [1] "gatito"

Ejercicio: Escribir una bucle que encuentre los cinemex alrededor de cada uno de los cinepolis de nuestro data.frame “cinepolis”, que convierta los data.frames resultantes de esa búsqueda en data.frames espaciales y que los guarde en disco con el nombre del cinepolis respectivo.

Cálculo de distancias entre puntos geográficos

R ofrece paquetes muy convenientes para calcular distancias y tiempos de recorrido entre localidades geográficas. Por ejemplo el paquete “gmapsdistance” permite calcular distancias y tiempos con distintos medios de transporte y partiendo en una fecha:hora determinada (sólo si se brina una key).

Por ejemplo podemos calcular la distancia y tiempo recorrido partiendo ahorita (es el default si no se especifica la fecha:hora de partida) al transportarse entre Toluca y CDMX:

library("gmapsdistance")
## Warning: package 'gmapsdistance' was built under R version 3.4.2
resultados = gmapsdistance(origin = "Toluca+Estado+De+Mexico+Mexico", 
                           destination = "CDMX+Mexico", 
                           mode = "driving")

resultados
## $Time
## [1] 4303
## 
## $Distance
## [1] 67749
## 
## $Status
## [1] "OK"

Esta función también entiende coordenadas lat lon, por ejemplo podemos encontrar la distancia y tiempo de recorrido caminando de la ciudad griega de Marathon a Atenas. Athens.

resultados = gmapsdistance(origin = "38.1621328+24.0029257",
                        destination = "37.9908372+23.7383394",
                        mode = "walking")

resultados
## $Time
## [1] 30487
## 
## $Distance
## [1] 39509
## 
## $Status
## [1] "OK"

Finalmente, se puede introducir a los argumentos “origin” y “destination” arreglos de localidades y se calcularán las combinaciones pertinentes.

or = c("Washington+DC", "New+York+NY", "Seattle+WA", "Miami+FL")
des = c("Los+Angeles+CA", "Austin+TX", "Chicago+IL", "Philadelphia+PA")

resultados = gmapsdistance(origin = or, destination = des, mode = "bicycling", combinations = "pairwise")

resultados
## $Time
##              or              de   Time
## 1 Washington+DC  Los+Angeles+CA 858830
## 2   New+York+NY       Austin+TX 598038
## 3    Seattle+WA      Chicago+IL 675002
## 4      Miami+FL Philadelphia+PA 411413
## 
## $Distance
##              or              de Distance
## 1 Washington+DC  Los+Angeles+CA  4582527
## 2   New+York+NY       Austin+TX  3138643
## 3    Seattle+WA      Chicago+IL  3587447
## 4      Miami+FL Philadelphia+PA  2173558
## 
## $Status
##              or              de status
## 1 Washington+DC  Los+Angeles+CA     OK
## 2   New+York+NY       Austin+TX     OK
## 3    Seattle+WA      Chicago+IL     OK
## 4      Miami+FL Philadelphia+PA     OK

Ejercicio: Calcular las distancias y tiempo de recorrido desde la CDMX a todos los puertos del país.

http://www.sct.gob.mx/puertos-y-marina/puertos-de-mexico/